NVU dynamics. I. Geodesic motion on the constant-potential-energy hypersurface 
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An algorithm is derived for computer simulation of geodesies on the constant potential-energy 
hypersurface of a system of N classical particles. First, a basic time-reversible geodesic algorithm 
is derived by discretizing the geodesic stationarity condition and implementing the constant poten- 
tial energy constraint via standard Lagrangian multipliers. The basic NVU algorithm is tested by 
single-precision computer simulations of the Lennard- Jones liquid. Excellent numerical stability is 
obtained if the force cutoff is smoothed and the two initial configurations have identical potential 
energy within machine precision. Nevertheless, just as for NVE algorithms, stabilizers are needed 
for very long runs in order to compensate for the accumulation of numerical errors that eventually 
lead to "entropic drift" of the potential energy towards higher values. A modification of the basic 
NVU algorithm is introduced that ensures potential-energy and step-length conservation; center- 
of-mass drift is also eliminated. Analytical arguments confirmed by simulations demonstrate that 
the modified NVU algorithm is absolutely stable. Finally, simulations show that the NVU algo- 
rithm and the standard leap-frog NVE algorithm have identical radial distribution functions for 
the Lennard- Jones liquid. 
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I. INTRODUCTION 

This paper and its companion Paper II [1| study NVU 
dynamics, i.e., dynamics that conserves the potential 
energy U for a system of N classical particles at con- 
stant volume V. NVU dynamics is deterministic and 
involves only the system's configurational degrees of free- 
dom. NVU dynamics is characterized by the system 
moving along a so-called geodesic curve on the constant 
potential-energy hypersurface Q defined by 



n = {(ri, r N ) G i? 3JV |f/(ri, r N ) = U } . (1) 

Mathematically, ft is a 3N — 1 dimensional differentiable 
manifold. Since it is imbedded in R 3N , SI has a natural 
Euclidean metric and it is thus a so-called Riemannian 
manifold Q . The differential geometry of hypersurfaces 
is discussed in, for instance, Ref. [jj 

A geodesic curve on a Riemannian manifold by defini- 
tion minimizes the distance between any two of its points 
that are sufficiently close to each other (the curve is char- 
acterized by realizing the "locally shortest distance" be- 
tween points). More generally, a geodesic is defined by 
the property that for any curve variation keeping the two 
end points Ha and fixed, to lowest order the curve 
length does not change, i.e., 
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Here dl denotes the line element of the metric. 

From a physical point of view it is sometimes useful 
to regard a geodesic as a curve along which the system 
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moves at constant velocity with zero friction. Such mo- 
tion means that at any time the force is perpendicular 
to the surface, and because the force performs no work, 
the kinetic energy is conserved. In this way geodesic 
motion generalizes Newton's first law, the law of iner- 
tia, to curved surfaces. The concept of geodesic mo- 
tion is central in general relativity, where motion in a 
gravitational field follows a geodesic curve in the four- 
dimensional curved space-time 

A general motivation for studying NVU dynamics is 
the following. Since all relevant information about a sys- 
tem is encoded in the potential-energy function, it is in- 
teresting from a philosophical point of view to study and 
compare different dynamics relating to U(ri, r^v). The 
"purest" of these dynamics does not involve momenta 
and relates only to configuration space. NVU dynam- 
ics provides such a dynamics. In contrast to Brownian 
dynamics, which also relates exclusively to the configu- 
rational degrees of freedom, NVU dynamics is determin- 
istic. NVU dynamics may be viewed as an attempt to 
understand the dynamic implications of the potential en- 
ergy landscape's geometryalong the lines of recent works 
by Stratt and coworkers 0, @] . 

Our interest in NVU dynamics originated in recent 
results concerning strongly correlating liquids and their 
isomorphs. A liquid is termed strongly correlating if 
there is more than 90% correlation between its virial and 
potential energy thermal equilibrium fluctuations in the 
NVT ensemble Q. The class of strongly correlating liq- 
uids includes most or all van der Waals and metallic liq- 
uids, whereas hydrogen-bonding, covalently bonded liq- 
uids, and ionic liquids are generally not strongly corre- 
lating. A liquid is strongly correlating if and only if it 
to a good approximation has "isomorphs" in its phase 
diagram p| Q • By definition two state points are isomor- 
phic [1[ if any two microconfigurations of the state points, 
which can be trivially scaled into one another, have iden- 
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tical canonical probabilities; an isomorph is a curve in the 
phase diagram for which any two pairs of points are iso- 
morphic. Only inverse-power-law liquids have exact iso- 
morphs, but simulations show that Lennard- Jones type 
liquids have isomorphs to a good approximation [§[ . This 
is consistent with these liquids being strongly correlating 
Q • Many properties are invariant along an isomorph, for 
instance the excess entropy, the isochoric heat capacity, 
scaled radial distribution functions, dynamic properties 
in reduced units, etc [1, the reduced- unit constant- 
potential-energy hypersurface Q is also invariant along 
an isomorph [8j . Given that several properties are invari- 
ant along a strongly correlating liquid's isomorphs and 
that H is invariant as well, an obvious idea is that Q's 
invariance is the fundamental fact from which all other 
isomorph invariants follow. For instance, the excess en- 
tropy is the logarithm of the area of Cl, so the excess 
entropy's isomorph invariance follows directly from that 
of 51. In order to understand the dynamic isomorph in- 
variants from the £1 perspective a dynamics is required 
that refers exclusively to Cl. One possibility is diffusive 
dynamics, but a mathematically even more elegant dy- 
namics on a diffcrcntiablc manifold is that of geodesies. 
Although these considerations were our original motiva- 
tion, it should be emphasized that that the concept of 
geodesic motion on 51 (or 51) is general and makes sense 
for any classical mechanical system, strongly correlating 
or not. 

We are not the first to consider dynamics on the 
constant-potential-energy hypersurface. In papers dating 
back to 1986 [l(| Cotterill and Madsen proposed a deter- 
ministic constant-potential-energy algorithm similar, but 
not identical, to the basic NVU algorithm derived below. 
Their algorithm was not discussed in relation to geodesic 
curves, but aimed at providing an alternative way to un- 
derstanding vacancy diffusion in crystals and, in partic- 
ular, to make easier the identification of energy barriers 
than from ordinary MD simulations. The latter prop- 
erty is not confirmed in the present papers, however 
we find that NVU dynamics in the thermodynamic limit 
becomes equivalent to standard NVE dynamics (Paper 
II [1]). Later Scala et al. studied diffusive dynamics on 
the constant-potential-energy hypersurface 51 [ll| , focus- 
ing on the entropic nature of barriers by regarding these 
as "bottlenecks". This point was also made by Cotter- 
ill and Madsen who viewed 51 as consisting of "pockets" 
connected by thin paths, referred to as "tubes", acting 
as entropy barriers. Reasoning along similar lines, Stratt 
and coworkers published in 2007 and 2010 three papers 
0, Q , which considered paths in the so-called potential- 
energy-landscape ensemble. This novel ensemble is de- 
fined as including all configurations with potential en- 
ergy less than or equal to some potential energy Uq. A 
geodesic in the potential-energy-landscape ensemble con- 
sists of a curve that is partly geodesic on the constant- 
potential-energy surface 51, partly a straight line in the 
space defined by U < Uq [5]. Stratt et aVs picture shifts 
"perspective from finding stationary points on the poten- 



tial energy landscape to finding and characterizing the 
accessible pathways through the landscape. Within this 
perspective pathways would be slow, not because they 
have to climb over high barriers, but because they have to 
take a long and tortuous route to avoid such barriers...." 
[5|. Thus the more "convoluted and laborinthine" the 
geodesies are, the slower is the dynamics Q. Apart from 
these three sources of inspiration to the present work, 
we note that geodesic motion on differentiable manifolds 
has been studied in several other contexts outside of pure 
mathematics, see, e.g., Ref. [l2|. 

The present paper derives and documents an algorithm 
for NVU geodesic dynamics. In Sec. II we derive the 
basic NVU algorithm. By construction this algorithm is 
time reversible, a feature that ensures a number of im- 
portant properties [l3|, [l4| • Section III discusses how to 
implement the NVU algorithm and tests improvements 
of the basic NVU algorithm designed for ensuring stabil- 
ity, which is done by single-precision simulations. This 
section arrives at the final NVU algorithm and demon- 
strates that it conserves potential energy, step length, 
and center-of-mass position in arbitrarily long simula- 
tions. Section IV briefly investigates the sampling prop- 
erties of the NVU algorithm, showing that it gives results 
for the Lennard- Jones liquid that are equivalent to those 
of standard NVE dynamics. Finally, Sec. V gives some 
concluding comments. Paper II compares NVU simula- 
tions to results for four other dynamics, concluding that 
NVU dynamics is a fully valid molecular dynamics. 



II. THE BASIC NVU ALGORITHM 

For simplicity of notation we consider in this paper 
only systems of particles of identical masses (Appendix A 
of Paper II generalizes the algorithm to systems of vary- 
ing particle masses) . The full set of positions in the 3N- 
dimensional configuration space is collectively denoted by 
R, i.e., 

R = (n,...,r w ). (3) 

Likewise, the full 3-/V-dimensional force vector is denoted 
by F. This section derives the basic NVU algorithm for 
geodesic motion on the constant-potential-energy hyper- 
surface 51 defined in Eq. ([T]) , an algorithm that allows one 
to compute the positions in step i + 1, Rj+i, from R,_i 
and R.;. Although a mathematical geodesic on a differ- 
entiable manifold is usually parameterized by its curve 
length [2], it is useful to think of a geodesic curve on 51 
as parameterized by time, and we shall refer to the steps 
of the algorithm as "time steps" . 

Locally, a geodesic is the shortest path between any 
two of its points. More precisely: 1) For any two points 
on a Riemannian manifold the shortest path between 
them is a geodesic; 2) The property of a curve being 
geodesic is locally defined; 3) a geodesic curve has the 
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property that for any two of its points, which are suf- 
ficiently close to each other, the curve gives the short- 
est path between them. A geodesic may, in fact, be the 
longest distance between two of its points. For instance, 
the shortest and the longest flight between two cities 
on our globe both follow great circles - these are both 
geodesies. In any case, the property of being geodesic is 
always equivalent to the curve length being stationary in 
the following sense: Small curve variations, which do not 
move the curve's end points, to lowest order lead to no 
change in the curve length. 

For motion on f2 the constraint of constant potential 
energy is taken into account by introducing Lagrangian 
multipliers. For each time step j there is the constraint 
U(Rj) = Uq and a corresponding Lagrangian multiplier 
Xj. Thus the stationarity condition Eq. @ for the dis- 
cretized curve length J2j |Rj ~ R-j— 1 1 subject to the con- 
straint of constant potential energy, is 



6 I £ |R, - Ri-xl - £ W R i)J = • ( 4 ) 

Since |R, • - R,-_i| = y/(Rj - R 3 -i) 2 and the 3N- 
dimensional force is given by Fj = —BU/dRj, putting 
to zero the variation of Eq. ((4]) with respect to Rj (i.e., 
the partial derivative d/dRi) leads to 

R-i — Ri— i Rz+i — Ri ^ ^ jp- o (5) 
|R» — R»-i| |R,+i — Rj| 

To solve these equations we make the ansatz of constant 
displacement length for each time step, 

|R,--R,-_i| =l (all j). (6) 

If the path discretization is thought of as defined by con- 
stant time increments, Eq. © corresponds to constant 
velocity in the geodesic motion. With this ansatz Eq. ([5]) 
becomes 

(Ri - Ri_i) + (Ri - Ri+i) + l \iFi = . (7) 

If &i = Hi — Hi i and bj = Rj — R,:+i, Eq. ([5]) implies 
a 2 = bf, i.e., = af-bf = (a^+bf) • (aj -bj). Since Eq. 
([7]) expresses that + b, is parallel to Fj , one concludes 
that Pi is perpendicular to a^ — b, = Rj+i — Rj-i- This 
implies 

Pi • Rj-i = Pi • Ri+l ■ (8) 

Taking the dot product of each side of Eq. ([7]) with Fj 
one gets 

F, • (R» - Ri-!) + F 4 • (Ri - R 1+1 ) + ZoAiF? = , (9) 



which via Eq. jSJ implies 

i oXi = _ 2 ^-(^_-Ri-0 , (10) 

i 

Substituting this into Eq. (JJJ) and isolating Rj+i we fi- 
nally arrive at 

R l+ i = 2R, - R,_i - 2[Fi • (Ri - R,_i)]F,/F 2 . (11) 

This equation determines a sequence of positions; it will 
be referred to as "the basic NVU algorithm" . The algo- 
rithm is initialized by choosing two nearby points in con- 
figuration space with the same potential energy within 
machine precision. 

The derivation of the basic NVU algorithm is com- 
pleted by checking its consistency with the constant 
step length ansatz Eq. ©: Rewriting Eq. (p~T|) as 
(Ri+i-Ri) = (R i -R i _i)-2[F r (R i - R,-_i)]Fj/F 2 we 
get by squaring each side (Rj+i — Ri) 2 = (Ri — Ri_i) 2 + 
4[F r (Ri - R l _ 1 )] 2 /Ff-4[F i -(R l - R,_i)] 2 /F 2 = (Rj- 
Ri_i) 2 . Thus the solution is consistent with the ansatz. 

Time reversibility of the basic NVU algorithm is 
checked by rewriting Eq. (fTT]) as follows 

Ri-i = 2Rj — Rj+i — 2[Fj • (Ri — Ri_i)]Fi/F 2 , (12) 
which via Eq. ((HJ implies 

Ri-i = 2Ri — Rj+i — 2[Pj • (Ri — Ri + i)]Fj/F 2 . (13) 

Comparing to Eq. (fTT]) shows that any se- 
quence of configurations generated by Eq. (|11[) 

Rj— i, Ri, Ri+i, ••• obeys Eq. (TTTj) in the time-reversed 
version R i+1 , R i; Ri_ l5 .... A more physical way to 
show that the basic NVU algorithm is time-reversal in- 
variant is to note that Eq. ([5]) is itself manifestly invari- 
ant if the indices i — 1 and i + 1 are interchanged. 

Appendix A shows that the basic NVU algorithm is 
symplectic, i.e., that it conserves the configuration-space 
volume element in the same way as NVE dynamics does. 
We finally consider potential-energy conservation in the 
basic NVU algorithm. A Taylor expansion implies via 
Eq. © that 

U l+1 - Ut-i = -F< ■ (Ri+i - Ri_i) + O(lg) - O(Zg) . 

(14) 

This ensures potential-energy conservation to a good ap- 
proximation if the discretization step is sufficiently small. 

The "potential energy contour tracing" (PECT) al- 
gorithm of Cotterill and Madsen [l(| is the following: 
Rj+i = 2Ri — Ri_i — [Fj • (Rj — Ri_i)]Fj/F 2 . Ex- 
cept for a factor of 2 this is identical to the basic NVU 
algorithm. The importance of this difference is appar- 
ent when it is realized that the PECT algorithm implies 
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Fj ■ (Rj+i — Rj) = 0, whereas it does not imply the time- 
reversed identity Fj • (Rj_i - R*) = 0. Thus the PECT 
algorithm is not time reversible. 

We end this section by reflecting on what is the relation 
between the NVU algorithm and continuous geodesic 
curves on Q. Can one expect that if the step length 
is decreased towards zero, the discrete sequence of points 
traced out by the algorithm converges to a continuous 
geodesic curve? The answer is yes, as is clear from the 
current literature that treats this problem in consider- 
able detail [15j. The literature deals with the analogous 
problem of classical mechanics where, as is well known, 
Newton's second law of motion can be derived from the 
principle of least action (Hamilton's principle). This is 
a variational principle. In the traditional approach one 
first derives continuous equations of motion from the vari- 
ational principle, then discretizes these equations to allow 
for computer simulations. Here we first discretized the 
quantity subject to the variational principle (Eq. (|4|) 
and only thereafter applied variational calculus. 

Euler himself first described discretization of time in 
the action integral, thus obtaining discretized versions 
of the Euler-Lagrange equations. There is now a large 
literature on this subject During the last decade, 

in particular, variational calculations applied after dis- 
cretization have come into focus in connection with for 
instance the development of algorithms for the control 
of robots. The general consensus is the following (we 
quote below from Ref. HH that provides an excellent, 
brief summary of the situation): "The driving idea be- 
hind this discrete geometric mechanics is to leverage the 
variational nature of mechanics and to preserve this vari- 
ational structure in the discrete setting... If one designs a 
discrete equivalent of the Lagrangian, then discrete equa- 
tions of motion can be easily derived from it by paral- 
leling the derivations followed in continuous case, and 
good numerical methods will come from discrete analogs 
to the Euler-Lagrange equations. In essence, good nu- 
merical methods will come from discrete analogs to the 
Euler-Lagrange equations - equations that truly derive 
from a variational principle... Results have been shown 
to be equal or superior to all other types of integrators 
for simulations of a large range of physical phenomenon, 
making this discrete geometric framework both versatile 
and powerful." 



III. TESTING AND IMPROVING THE BASIC 
NVU ALGORITHM 



This section discusses the numerical implementation of 
the basic NVU algorithm and how to deal with round- 
off errors that arise for very long simulations. The model 
system studied is the standard Lennard- Jones (LJ) liq- 
uid with N — 1024 particles. Recall that the LJ pair 
potential v(r) is given by 



v(r) 



As 







r 



(15) 



Here e sets the energy scale and a the length scale; hence- 
forth the unit system is adopted in which these quantities 
are both unity. All simulations except those of Fig. pre- 
fer to the state point with density 0.85 and temperature 
0.7 in reduced units. The two initial configurations were 
taken from NVE simulations of this state point. Unless 
otherwise specified the forces and their derivative were 
adjusted to be continuous via smoothing from a value 
just below the cutoff distance r c to r c . We refer to this 
as a "smoothed force potential" . The cutoff distance was 
chosen as the standard LJ cutoff r c = 2.5o\ The simula- 
tions were performed using periodic boundary conditions. 
In order to easier test the numerical stability of the NVU 
algorithm, simulations were performed in single precision 

m. 



A. Implementing the basic NVU algorithm 

We rewrite Eq. (|11[) into a leap-frog version by intro- 
ducing new variables defined by 



Ai+i/a = Rm - R* . (16) 
In terms of these variables the basic NVU algorithm is 

A l+1/2 = Ai_ 1/2 - 2(F< • Ai_ 1/2 )F i /F' 
Rj+i = R s : + Aj +1 / 2 . (17) 

The equations (fTTj) are formally equivalent to Eq. (jlll) . 
Numerically, however, they are not equivalent and - as 
is also the case for standard NVE dynamics - the leap- 
frog version is preferable because it deals with position 
changes [iH ]. 

Figure[T](a) shows the potential energy as a function of 
time-step number. The system's potential energy jumps 
every second step, jumping between two distinct values 
(inset). This is also reflected in the distribution of the 
quantity IqXi shown in green in Fig. [3c). A priori one 
would expect a Gaussian single-peak distribution of loXi, 
but the distribution has two peaks. What causes the 
potential energy to zig-zag in an algorithm constructed 
to conserve potential energy? The answer is hinted at 
in Eq. (TTJJ according to which the NVU algorithm im- 
plies energy conservation to a good accuracy, but only 
every second step. Thus if the two initial configurations 
do not have identical potential energy, the potential en- 
ergy will zig-zag between two values. Figure QJb) shows 
that even if a simulation is initiated from two configura- 
tions with very close potential energies, the zig-zag phe- 
nomenon persists, though now on a much smaller scale. 




4.0x10"' 6.0X10" 6 8.0X10" 6 l.OxlO" 5 1.2x10" 



FIG. 1: (a) Evolution of the potential energy U according 
to the basic NVU algorithm (Eq. (|17[) ) started from two 
consecutive configurations of an NVE simulation. The inset 
shows a snapshot of the first ten integration steps where lines 
connect the data points; clearly the system jumps distinctly 
between two potential-energy hypersurfaces. (b) Evolution 
of U started from two configurations with a very small po- 
tential energy difference. The algorithm still jumps between 
two potential-energy hypersurfaces, but the difference is much 
smaller, (c) Probability distribution of the Lagrangian mul- 
tiplier times the length lo, lo\ of Eq. (|I0[) . obtained from 
simulations over 2.5 TO 6 steps. The green distribution corre- 
sponds to (a), the blue distribution to (b). 



There are further numerical issues that effect the sta- 
bility of the basic NVU algorithm. In Fig. [5] the evo- 
lution of the potential energy is given for a long simula- 
tion, which also includes data from simulations using a 
smoothed force potential. Better numerical stability is 
clearly obtained for the smoothed force potential (black 
curve), but smoothing does not ensure a constant poten- 
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FIG. 2: Evolution of \(U -U(0))/U(0)\ for a simulation using 
the basic NVU algorithm. The red curve gives results from a 
simulation where the potential is cut and shifted at r = 2.5a, 
the black curve gives results for a smoothed force potential. 

tial energy and absolute stability. 



B. Improving the algorithm to conserve potential 
energy and step length indefinitely 

The last subsection showed that using a smoothed force 
potential and ensuring that the two starting configura- 
tions have identical potential energy within machine pre- 
cision, a more stable algorithm is arrived at. Neverthe- 
less, absolute stability is not obtained. This is illustrated 
in Fig. Eta), which shows that the potential energy for 
a system with a smoothed force potential over five mil- 
lion time steps still exhibits a slight "entropic drift" (red 
curve) . By entropic drift we mean the drift due to round- 
off errors - a drift that unavoidably takes the system to 
higher energies because there are many more such states, 
an entropic effect. Figure [3Jb) shows that also the step 
length is not conserved. Both problems are caused by the 
accumulation of round-off errors. These problems are less 
severe if one switches to double precision, of course, but 
for long simulations entropic drift eventually sets in (for 
billions of time steps). 

We would like to have an algorithm that is absolutely 
stable, i.e., one that does not allow for any long-time drift 
of quantities the basic NVU algorithm was constructed 
to conserve: the potential energy, the step length, and 
the center of mass (CM) position (just as in standard 
NVE dynamics the CM position is exactly conserved in 
the basic NVU algorithm Eq. (jTTJ) because the forces 
sum to zero due to the translational invariance of the 
potential energy: U(ri +r°, r« + r°) = J7(ri , rjv))- 

Drift of the CM position is trivially eliminated by ad- 
justing the particle displacements according to Ar„ = 
Ar„ — J2 n Ar n /N, every 100'th time steps. This correc- 
tion corresponds to setting to zero the total momentum 
of the system in an NVE simulation. 

Robust potential energy conservation is obtained by 
adding a term that is zero if the potential energy equals 
the target potential energy U (this quantity was previ- 
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ously denoted by Uq, but to avoid confusion with the 
time step index we drop the subscript zero), 



-1/2 



-1/2- 



- 2F, ; • A, 



-1/2- 



-c/.-i-ur./F 



(18) 

To show that this modification of the NVU algorithm 
prevents drift of the potential energy, we take the dot 
product of each side of Eq. ([T5]) with Fj, leading to 



L i+l/2 



Ai_i /2 ) = Ui-x - 
Fj • (Rj+i — Rj_i) 



•A 

- U. 

= -(U i+1 -Ui 



-1/2 + Ui. 

Since F, 



i-Uor Fi-(A i+1/2 + 
+ Ai_ 1/2 ) = 
), this implies 



-i) + O(l 3 



Ui. 



\J + 0(l 3 Q ) 



(19) 



Thus entropic drift has been eliminated and the potential 
energy is conserved indefinitely except for small fluctua- 
tions. 

We next address the problem of conserving step length. 
This is ensured by the following modification of the al- 
gorithm, 



-1/2 



lo 



Ai_i /2 + (-2P 4 • A,_ 1/2 + - U)Fi/F? 



-1/2 



-2F< ■ A,- 



-1/2 



U-! 



U)Pi/Ff| 
(20) 

Equation ([20]) gives what we term the final NVU algo- 
rithm (for brevity: "the NVU algorithm", in contrast 
to Eq. (fTTj) that is referred to as "the basic NVU algo- 
rithm" ) . 

In simulations the NVU algorithm is implemented as 
follows. The target potential energy U is chosen from 
an NVE or an NVT simulation at the relevant state 
point. The step length Iq is chosen according to the ac- 
curacy aimed for. Suppose at a given time the quantities 
Hi, Aj_!/ 2 , and Ui-! are given. From Rj the forces Fj 

Fj, and Ui—! the quantity 
(l20l) . Finally, the positions 



1/2; 



are calculated. From A 
A i+1 / 2 is calculated via Eq 
are updated via H4+1 = R 



Ai+1/2 an d the potential 



energy is updated via U — E/(Rj). 

By construction the NVU algorithm Eq. (l2Tfl) ensures 
constant step length, 



l A i+l/2 



lo ■ 



(21) 



but is the potential energy still conserved for arbi- 
trarily long runs? If the denominator of Eq. ([20]) 
is denoted by Di, taking the dot product of each 
side of this equation with Fj leads to Fj • A i+1 / 2 = 
(Zo/A) [-Fj • Aj_ 1/2 + U-! - U] . Writing h/D t = 1 + 
Si in which Si — 0(1$) with p > 1, we get Fj • (A i+1 / 2 



^i-1/2) 
since Fj 
Fj-Aj 



S. t [-Fj ■ Ai- 1/2 ] + (1 + Si) [U-! - U]. Thus, 



-1/2 



(Aj + i /2 + Aj_ 1/2 ) - U-! - Ui +1 + 0(l 3 ) and 
= Ui-!-Ui+0(l%), we get V-Ui+i+Oitf) = 



Si [Ui — U + O(Zq)] . This implies again 



U i+ ! = U + 0(ig) 



(22) 
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FIG. 3: (a) Evolution of U with and without the numerical 
stabilization (Eq. (|20p ~): The red curve gives results using 
the basic NVU algorithm Eq. ()17p with two identical initial 
potential energies and smoothed force potential. The black 
curve gives simulation results under the same conditions using 
the final NVU algorithm (Eq. (20])). (b) Evolution of the step 
length I Ai| for the same simulations. 



In summary, for simulations of indefinite length the 
NVU algorithm Eq. (l2Tfl) ensures constant step length 
and avoids entropic drift of the potential energy. Figure 
[3Ja) shows the evolution of the potential energy using 
the basic NVU algorithm (red) and the final NVU algo- 
rithm (black), Fig. |3Jb) shows the analogous step length 
evolution. Figure SJa) shows that the distribution of the 
Lagrangian multiplier is only slightly affected by going 
from the basic (red) to the final (black) NVU algorithm. 
Figure IHh) shows the evolution of Si in the final NVU 
algorithm, which as expected is close to zero. 

We remind the reader that the modifications were in- 
troduced to compensate for the effects of accumulating 
random numerical errors for very long runs, and that the 
modifications introduced in the final NVU algorithm Eq. 
(|20p vanish numerically in the mean. The prize paid for 
stabilizing the basic NVU algorithm is that the full NVU 
algorithm is not time reversible. In view of the fact that 
the improvements introduced to ensure stability lead to 
very small corrections, the (regrettable) fact that the cor- 
rections violate time reversibility is not important. 
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" s lu 200 400 600 800 1000 
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FIG. 4: (a) The distribution of the Lagrangian multiplier 
times lo with (black) and without (red) the numerical stabi- 
lization of the final NVU algorithm Eq. p0|. (b) Evolution 
of the quantity 5, defined by lo/Di = I + Si ; as expected this 
quantity is small and averages to zero. 



IV. SAMPLING PROPERTIES OF THE NVU 
ALGORITHM 

In order to investigate whether the NVU algorithm 
gives physically reasonable results we compare results 
from NVU and NVE simulations for the average of a 
quantity that depends only on configurational degrees of 
freedom. This is done in Fig. [5j which shows the radial 
distribution function g(r) at three state points. The red 
dots give NVU simulation results, the black curve NVE. 
Clearly the two algorithms give the same results. This 
finding is consistent with the conjecture that the NVU 
algorithm probes all points on Q with equal probabil- 
ity. Note that this is not mathematically equivalent to 
conjecturing that the NVU algorithm probes the config- 
uration space microcanonical ensemble, which has equal 
probability density everywhere in a thin energy shell be- 
tween a pair of constant-potential-energy manifolds. The 
latter distribution would imply a density of points on Q 
proportional to the length of the gradient of U(R) (the 
force) , but this distribution cannot be the correct equilib- 
rium distribution because the basic NVU algorithm Eq. 
(fTTj) is invariant to local scaling of the force. In the ther- 
modynamic limit, however, the length of the force vector 
becomes almost constant and the difference between the 
configuration space microcanonical ensemble and the f2 



s- 
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FIG. 5: Radial distribution functions g(r) for a single com- 
ponent Lennard- Jones system at the following state points: 
(a) T = 2.32 and p = 0.85; (b) T = 1.1 and p = 0.427; 
(c) the crystal at T — 0.28 and p — 0.85. The black curves 
show results from NVE simulations, the red dots from NVU 
simulations (Eq. ([20}). 

equal- measure ensemble becomes insignificant. 

Paper II details a comparison of NVU dynamics to 
four other dynamics, including two stochastic dynam- 
ics. Here both simulation and theory lead to the con- 
clusion that NVU and NVE dynamics are equivalent in 
the thermodynamic limit. 

V. CONCLUDING REMARKS 

An algorithm for geodesic motion on the constant- 
potential-energy hypersurface has been developed (Eq. 
(P2P1) ). Single-precision simulations show that this algo- 



8 



rithm, in conjunction with compensation for center-of- 
mass drift, is absolutely stable in the sense that poten- 
tial energy, step length, and center-of-mass position are 
conserved for indefinitely long runs. The algorithm re- 
produces the NVE radial distribution function of the 
LJ liquid, strongly indicating that correct configuration- 
space averages are arrived at in NVU dynamics. 

Although NVU dynamics has no kinetic energy pro- 
viding a heat bath, it does allow for a realistic description 
of processes that are unlikely because they are thermally 
activated with energy barriers that are large compared to 
fcgT (Paper II). In NVU dynamics, whenever a molec- 
ular rearrangement requires excess energy to accumulate 
locally, this extra energy is provided by the surrounding 
configurational degrees of freedom. These provide a heat 
bath in much the same way as the kinetic energy provides 
a heat bath for standard Newtonian NVE dynamics. 

The companion Paper (II) compares the dynamics 
of the Kob-Andersen binary Lennard- Jones liquid sim- 
ulated by the NVU algorithm and four other algorithms 
(NVE, NVT, diffusion on Q, Monte Carlo dynamics), 
concluding that results are equivalent for the slow de- 
grees of freedom. Paper II further argues from simula- 
tions and non-rigorous argumens that NVU dynamics 
becomes equivalent to NVE dynamics as N — > oo. 



The Jacobian of this map J(Si — > S^+i) is given by 



2-2 



FiCEti-Ri-i) 



-2- 



; , K " p ," / 

F i (H. i -¥t. i _ 1 ) 



-1 + 2- 



1 

o 



2-2- 



9X2, i 



(A3) 

This may be regarded as a two-by-two block matrix con- 
siting of blocks A,B,C,D. The determinant of this 
block matrix is |J| = |AD-BC| = |-BC| = (-1) M |B|, 
giving (where the index i is dropped for brevity) 
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Appendix A: Proof that the basic NVU algorithm is 
symplectic 

This Appendix proves that the basic NVU algorithm 
conserves the configuration-space volume element on the 
hypersurface O in the same sense as the NVE algorithm 
conserves the configuration-space volume element. We 
view the basic NVU algorithm (Eq. (fTTjl h 



Ri+i — 2R; — Rj_i 2 -F, , (Al) 



|J| = (-1) 
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-1 + 2^ 
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(A4) 

Defining the unit- length vector n = /F 2 , the last equality 
of Eq. (TA4| follows from B = 1 + 2n • n T => B 2 = 
l + 4n n T -4n n T = 1. Since \B\ 2 = \B 2 \ = 1, one has 
\B\ = ±1. Thus the volume element transforms as 



as a mapping of R 6N into itself. In the 6N dimensional 
configuration space of subsequent time-step pairs Si = 
{Ri,R,_i}, the NVU algorithm is 



(A5) 



Si — > Sj+l — {R;+i, R;} — {2Ri — R,_l 



2Fi • (Rj — Rj_ 



This means that the basic NVU algorithm conserves 
the volume element in the 6N dimensional configuration 
space, i.e., that the algorithm is symplectic just as the 
i4^.Ra)&ai)ithm is. 
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